---
title: "Dead"
output: github_document
---

```{r loading, eval = T}

library(readr)
library(dplyr)
library(tidyr)
library(zoo)
library(data.table)
library(parallel)
library(gridExtra)
library(ggplot2)
library(knitr)

raw_dat.df <- 
  read_csv("istat-5-mar-2021/comuni_giornaliero_31dicembre.csv.zip", 
           locale = locale(encoding = "ISO-8859-1")) %>%
  select(!(REG:NOME_COMUNE) & !(TIPO_COMUNE) & !(M_11:F_20))
  
geo_ids.v <- 
  unique(raw_dat.df$COD_PROVCOM)

to_include <- unique(raw_dat.df$COD_PROVCOM[raw_dat.df$T_20 != "n.d."])

length(to_include) == length(geo_ids.v)

all_days <- unique(raw_dat.df$GE)
all_days <- all_days[grepl("^01|^02|^03|^04|^05|^06|^07|^08|^09|^10|^11|^12", all_days)]
all_days <- all_days[order(all_days)]

```

```{r slicing, eval = F}

# Slice 

splitVec <- 
  function(x, n) split(x, cut(seq_along(x), n, labels = FALSE)) 

geo_ids.list <- 
  splitVec(geo_ids.v, 100)

dat_over65.list <- list()

for (i in 1:length(geo_ids.list)) {
  
  print(i)
  
  m <- 
    expand.grid(geo_ids.list[[i]], 
                all_days, 
                paste0("T_", 15:20), 
                unique(raw_dat.df$CL_ETA))
  colnames(m) <- c("COD_PROVCOM","GE", "year", "CL_ETA")
  
  dt <- 
    data.table(m, key=colnames(m))
  
  dat.dt <- 
    raw_dat.df %>%
    dplyr::filter(CL_ETA >= 14) %>% # Over 65
    dplyr::filter(COD_PROVCOM %in% geo_ids.list[[i]]) %>%
    tidyr::gather(year, value, T_15:T_20) %>%
    data.table(key=colnames(m))
  
  rm(m)
  
  dat.dt <-
    merge(dat.dt, dt, all = T)
  dat.dt$value[is.na(dat.dt$value)] <- "0"
  
  rm(dt)
  
  dat_over65.dt <- 
    dat.dt %>%
    dplyr::mutate(value = as.numeric(value)) %>%
    dplyr::group_by(COD_PROVCOM, GE, year) %>%
    dplyr::summarize(value = sum(value, na.rm = T), .groups = 'drop')
  
  rm(dat.dt)
  
  dat_over65.list[[i]] <-
    dat_over65.dt %>%
    dplyr::group_by(COD_PROVCOM, year) %>%
    dplyr::arrange(as.character(GE)) %>%
    dplyr::mutate(value_ma14 = rollmean(value, 14, align = 'right', fill = NA))
  
  rm(dat_over65.dt)
  
}

```

```{r rbindlist, eval=F}

dat_over65.dt <- 
  rbindlist(dat_over65.list)

rm(dat_over65.list)

save(all_days, dat_over65.dt, file = "dat_over65.dt.RData")

```


```{r eval=T}

library(sf)

load("dat_over65.dt.RData")

load("../grid/hex_ita_10km.sf.RData")
st_crs(hex_ita_10km.sf) <- st_crs(32632)

# hex_ita_1km_union.sf <- 
#   st_union(hex_ita_10km.sf %>% dplyr::filter(!is.na(CODREG)))

italy.sf <- 
  read_sf("../shapefile/Italy/Limiti01012020/Com01012020/Com01012020_clipped_to_hex.shp") %>% 
  sf::st_simplify(preserveTopology = TRUE, dTolerance = 1000) %>%
  sf::st_make_valid()

ggplot(italy.sf) + 
  geom_sf(aes(fill = PRO_COM_T %in% to_include), colour = NA) +
  scale_fill_grey() + theme_bw() + labs(fill = "Data available")

raw_dat.df %>%
  dplyr::summarize(`pop. covered %` = sum(T_19[CL_ETA >= 15 & 
                       COD_PROVCOM %in% to_include]  / 
                  sum(T_19[CL_ETA >= 15]))) %>%
  kable()
  
# italy.sf <-
#   italy.sf %>%
#   sf::st_simplify(preserveTopology = TRUE, dTolerance = 1000) %>%
#   st_make_valid()

dat_over65_spread.dt <- 
  dat_over65.dt %>%
  dplyr::select(-value_ma14) %>%
  tidyr::pivot_wider(names_from = year, values_from = value)

dat_over65_spread.dt$GE <- 
  as.character(dat_over65_spread.dt$GE)
dat_over65_spread.dt$COD_PROVCOM <-
  as.character(dat_over65_spread.dt$COD_PROVCOM)

fun1 <- function(day) {
  this.dt <- 
    dat_over65_spread.dt %>% 
    dplyr::filter(GE == day)
  this.sf <- 
    merge(italy.sf, this.dt, 
          by.x = "PRO_COM_T", by.y = "COD_PROVCOM")
  return(this.sf)
}

these.sf <- lapply(all_days, fun1)
names(these.sf) <- all_days

fun2 <- function(i, these.sf, hex_ita_10km.sf) {
  this.sf <- these.sf[[i]]
  require(sf)
  hex_dead.df <- data.frame()
  for (year in c(paste0("T_", 15:20))) {
      this_sampling.sf <- 
        this.sf[this.sf[[year]]>0,]
      if (nrow(this_sampling.sf) == 0) {
        next
      }
      this_sample.sf <- 
        st_sample(this_sampling.sf, this_sampling.sf[[year]])
      this_res <- 
        as.data.frame(table(unlist(st_within(this_sample.sf, hex_ita_10km.sf))))
      
      this_res$year <- year
      this_res$day <- names(these.sf)[[i]]
      hex_dead.df <- rbind(hex_dead.df, this_res)
  }
  return(hex_dead.df)
}

cl <- makeCluster(10)

hex_dead.list <- parLapply(cl, 1:length(all_days), fun2, these.sf, hex_ita_10km.sf)

stopCluster(cl)

hex_dead.df <- bind_rows(hex_dead.list)
colnames(hex_dead.df)[1] <- "hex_id"
save(hex_dead.df, file = "hex_dead.df.RData")
```

```{r}
hex_dead.df %>%
  dplyr::group_by(day)
```


```{r, fig.width = 9, fig.height=13}

plotFun <- function(day) {
  this.sf <- hex_ita_10km.sf
  this.sf$value <-
    hex_dead.df[hex_dead.df$day == day,]$Freq[match(this.sf$hex_id,
                                                    hex_dead.df[hex_dead.df$day == day,]$hex_id)]
  this.sf$value[is.na(this.sf$value)] <- 0
  return(ggplot(this.sf) + 
           geom_sf(aes(fill=value), colour = NA) + 
           scale_fill_distiller(palette="Spectral", direction = -1) +
           labs(fill = "Dead", caption = day) +
           theme(axis.text = element_blank(),
                 axis.ticks = element_blank()))
}

plots <- lapply(c("0101", "0115", "0131",
                  "0201", "0215", "0229",
                  "0301", "0315", "0331",
                  "1001", "1101", "1201"), 
                plotFun)



do.call("grid.arrange", c(plots, ncol=3))

```
